Second order tunneling of two interacting bosons in a driven triple well 
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We investigate quantum tunneling of two repulsive bosons in a triple-well potential subject to a 
high-frequency driving field. By means of the multiple-time-scale asymptotic analysis, we evidence 
a far-resonant strongly-interacting regime in which the selected coherent destruction of tunneling 
can occur between the paired states and unpaired states, and the dominant tunneling of the paired 
states is a second order process. Two Floquet quasienergy bands of the both kinds of states are given 
analytically, where a fine structure up to the second order corrections is displayed. The analytical 
results are confirmed numerically based on the exact model, and may be particularly relevant to 
controlling correlated tunneling in experiments. 
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I. INTRODUCTION 

Advances in laser technology have enabled studies of 
quantum tunneling and its coherent control for a single 
particle in light-induced quantum wells without dissipa- 
tion jlj]. Research attempting to manipulate quantum 
states has been underway for a long time [2, S|- The 
time-periodic driving field is a powerful tool to control 
the tunneling dynamics and can lead to important phe- 
nomena, such as dynamic localization (DL) [J, [a], coher- 
ent destruction of tunneling (CDT) [1,01: and photon- 
assisted tunneling |8l-[l0|. In recent years, the effects of 
interparticle interaction have attracted much attention. 
It was shown that adjusting the interaction can give rise 
to richer behavior, including many-body selective CDT 
prl [13 ] and the second order tunneling of two interacting 
bosons [ij|. The two-body interaction model is the sim- 
plest model for studying the interacting effects, and has 
received much attention Il4l4l7l , since the seminal exper- 
imental result was reported [18| . The tunneling dynamics 
is related to the interplay between the interparticle inter- 
action and external field, and the former can be tuned by 
the Feshbach resonance technique p^ . 

In the presence of interaction and periodic external 
field, the quantum well system may be nonintegrable that 
necessitates the perturbation method for an analytical 
investigation. The multiple-scale technique is a very use- 
ful perturbation method and has been extensively em- 
ployed for different physical systems |20l - |24| . It was 
demonstrated that with the multiple-scale perturbation 
method, the usual high-frequency approximation corre- 
sponds to the first-order perturbation correction [2J]. 
In the far-resonant strongly-interacting regime with a 
stronger reduced interaction |25j , the high-frequency ap- 
proximation is no longer valid. In this case, the dominant 
tunneling of paired states is a second-order process of 
long time scale and it can be described by the second- 
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order perturbation correction. The correlated tunnel- 
ing of two strongly-interacting atoms corresponding to 
time-resolved second order tunneling has been observed 
directly in an undriven double well system [l3|. Very 
recently, Longhi et al studied the second order effect of 
two far-resonant strongly-interacting bosons in a periodi- 
cally driven optical lattice by using a multiple-time-scale 
asymptotic analysis [26J. 

As above-mentioned, a lot of works on tunneling dy- 
namics of two interacting atoms focus on the systems 
with double-well or optical lattice potentials. The triple- 
well system is a bridge between the double well and the 
optical lattice systems, and is very important for us to 
fully understand coherent control of particle tunneling 
in the quantum wells |27h33 |. Besides, the triple- well 
system itself owns some novel phenomena, e.g., the stim- 
ulated Raman adiabatic passage ;27], which is a scheme 
that adiabatically transport a quantum particle from the 
left well to the right well with negligible middle well oc- 
cupation at all times. The tunneling dynamics of two- 
particle in a triple- well system have also attracted exten- 
sive attention ^33;h35i] . however, research on the second 
order effect of the system has not been reported yet. 

In this paper, we investigate the coherent control of the 
second order tunneling for two triple- well confined bosons 
driven by a high-frequency laser field. By means of the 
multiple-time-scale asymptotic analysis, we characterize 
quantum dynamics of the two bosons with the contin- 
uous increase of interaction intensity, and demonstrate 
a far-resonant strongly-interacting regime in which two 
bosons initially occupying the same well would form a 
stable bound pair, because of the selected CDT between 
the paired states and the unpaired states. Taking into ac- 
count the second order tunneling effect, the prediction on 
the CDT is confirmed by the Floquet quasienergy anal- 
ysis, where the Floquet quasienergy band of the three 
unpaired states exhibits the avoided level-crossings (or 
new level-crossings) at (or near) the collapse points, and 
the fine structure of quasienergy band of the three paired 
states shows the different level-crossings beyond the for- 
mer collapse points. Good agreements between the ana- 



lytical and numerical results are shown, which could be 
verified further under the current accessible experimental 
setups [l3,[ll[3il. 



II. THE MODEL AND HIGH-FREQUENCY 
APPROXIMATION 



We consider two interacting bosons confined in a triple- 
well potential and driven by an ac field. The Hamiltonian 
of the system in the tight-binding approximation is de- 
scribed by the three-site Bose-Hubbard model |33l435| 



H{t) = - J{a\a2 
. Uo 
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where the operator aj annihilates (creates) a boson 
in well I; J denotes the nearest- neighbor hopping ma- 
trix element, Uq is the on-site interaction energy, and 
e{t) = ecos(cjt) is the ac driving of amplitude e and fre- 
quency UJ. For simplicity, we adopt h = \ throughout this 
paper. The reference frequency wq = lOOHz is used to 
normalize the energy and the parameters J, Uq, s and w. 



such that all the 
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and time t is normalized in units of w, 
quantities become dimensionless 37|, 
assumed that the three wells are deep enough such that 
the Wannier functions of the two interacting bosons be- 
longing to different wells have very small overlap. A Fock 
basis \Nl, Nm^Nh) is useful to describe the two interact- 
ing bosons in the triple- well system, where N^, Nm and 
Nn are the number of bosons localized in the left, middle 
and right wells, respectively, with Nl + Nm + Nn = 2. 
The quantum state |V'(0) of ^^^ system is expanded as 
the linear superposition of the Fock states, 

m)) = ci(i)|2, 0, 0) + C2(t)|0, 2, 0) + C3(i)|0, 0, 2) + 

C4(i)|l, 1,0) + C5(t)|l,0, 1) + ceit)\0, 1, 1), (2) 

where Cj{t) [j = 1,2,..., 6) denote the time-dependent 
probability amplitudes of finding the two bosons in the 
six different Fock states and they obey the normaliza- 
tion condition X]7=i ki(OP = 1- Inserting Eqs. (1) and 
(2) into Schrodinger equation idt\4'{t)) — H{t)\ip{t)), one 
obtains the coupled equations for the amplitudes Cj{t) 

ici ^ [Uo - 2e(i)]ci - \/2Jc4, 

ic2 = U0C2 - \/2J(c4 + ce), 

ic3^[Uo + 2e{t)]c3-V2Jce, 

ici = -e{t)ci - J(V2ci -I- V2c2 + C5), 

ic5 = ~J{ci + ce), 

ice = e{t)ce - J{\f2c2 + \/2c3 + C5). (3) 

Although it is difficult to obtain exact analytical solu- 
tions of Eq. (3), we can approximately study some in- 
teresting phenomena in the high-frequency regime with 



a; ^ J. To do so, we rewrite the interaction strength as 
t/o — 'Tiw -I- u for |w| < w/2, 771 = 0, 1, 2, ... with u being 
the reduced interaction strength [23 , and make the func- 
tion transformations ci(t) = ai(t) exp[— zC/q^ -I- 2i(/?(t)], 
02(1) = a2(t)exp(-iC/o^), C3(i) = a3(t) exp[-iC/ot - 
2iiy9(i)], C4(i) — a4(t)exp[i(^(i)], c^ii) — a^{i), and 
C6(i) = a6(i) exp[— iiy9(i)], with aj(t) being the slowly- 
varying functions and ip(i) — L s cos{uJT)dT = — sin(a;t). 
Then, Eq. (3) is transformed into the coupled equations 
in terms of aj{t). Under the high-frequency approxi- 
mation, the rapidly oscillating functions included in the 
equations can be replaced by their time average such that 
the equations of Qj (t) become j35| 
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where J7m is the TTith-order Bessel function of the first 
kind, and e**"* are the slowly varying functions for a 
small u value. 

In Ref. [2^, Longhi proposed that the well-known 
high-frequency approximation commonly used to study 
CDT corresponds to the first-order perturbation approx- 
imation of the multiple-time-scale asymptotic analysis. 
If the first-order correction term vanishes in the pertur- 
bation treatment, the high-order corrected terms become 
important. Noticing that for a set of fixed external field 
parameters dynamical behavior of the system (j4]) is re- 
lated to the self-interaction intensity. In this work, we 
do not concern about the very strong interaction (e.g., 
Uo > 6a;), since for such a interaction we need to consider 
not only the usual on-site atom-interaction strength, but 
also the interactions between atoms on neighboring lat- 
tice sites [33, which is beyond the considered case. 

When the condition J'o = is satisfied, Eq. (|4]) 
shows that CDT occurs for the weakly-interacting case 
{m = 0, \u\ <C a;) that leads all the first derivatives of the 
probability amplitudes to zero. This can be further con- 
firmed by calculation of the Floquet quasienergies of the 
system in Sec. IV. Besides, for the resonant strongly- 
interacting case (w = 0,TO = 1,2,...), CDT for paired 
states is observed when the condition Jm = is satis- 
fied that leads the first derivatives aj{t) (j — 1,2,3) of 
paired-state amplitudes to zero. This is consistent with 
that of two interacting electrons in quantum dot arrays 
by numerical computation of the Floquet quasienergies 



|40| . As an example, we show time evolutions of the 
probabilities Pj(t) = jc^p — |ajP(j ~ 1,2,..., 6) for the 
resonant case with J = 1, Uq =^ cu ^ 80, e/cu — 2.405 
and P2(0) = 1, Pj=i2{0) — 0, as in Fig. 1, where the first 
order result (the circular points) from Eq. (4) is con- 
firmed by the direct numerical simulation (the curves) of 
Eq. (3). From Fig. 1, we can see that transitions be- 
tween the paired states with probabilities P,-, j = 1,2,3 
in Fig. 1(a) and the unpaired states with probabilities 
Pj, j = 4,6 in Fig. 1(b) happen periodically for the res- 
onant case. We will come back to this property for com- 
parison with the difference from the far-resonant case in 
next section. 




FIG. 1: (Color online) Time evolutions of the probabilities Pj — 
l^il (J — 1' 2, ..., 6) in six different states for two bosons initially occu- 
pying the middle well. The parameters are set as J — 1, Uq — uj — SO, 
e/u} — 2.405. (a) The probabilities of the three paired states, where the 
dashed line corresponds to P2 , and the solid line to Pi, 3. (b) The prob- 
abilities of the three unpaired states, where the dashed line associates 
with P5, and the solid line with P4,6. The circular points indicate the 
numerical results from the first approximate Eq. (4) and the curves 
describe the numerical solutions of the original Eq. (3). Hereafter, all 
variables and parameters are dimensionless. 



time. The probability amplitudes aj{t') {j — 1,2, ...,6) 
are expanded as a power series of e 
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Owing to the high order infinitesimal can be neglected 
in the high-frequency regime, we approximately rewrite 
the probability amplitudes as the leading order aj{t') — 
af\t') = A,it'). Thus, |^,|2 = |c,f (j = 1,2,...,6) 
denote the probabilities of finding the two bosons in the 
six different Fock states in Eq. (2). According to the 
perturbation analysis in the Appendix, we readily obtain 
that such amplitudes are the slowly- varying functions in 
time, which satisfy the following linear equations with 
constant coefficients, 

■ dA2 ^2 



^^ = 2e^ [2^2/91 + (Ai + ^3)^2], 

dA^ „ 9, , , . 

i—f ^ 2e^A3Pi + A2P2); 



(6) 



dt' 
«^ = -eJo(^)A5 - 2e^2Aipi + ^6^2), 

^^ = -eM-){A,+Ae), 
dt' uj 

i^ = -eJo(-)A5 - 2e\2A,p^ + A4P2), (7) 
dt' u) 

where pi ( i = 1, 2) are set as (see Appendix) 



We have known that Eq. (4) is a good approximation 
of Eq. (3) only for small values of the reduced interaction 
strength, |u| ^ lo. When the |u| values tend to their max- 
imum |m| — a;/2, the functions e*™* vary middlingly fast 
compared to the rapidly oscillating driving field. Con- 
sequently, in Eq. (4), although e^™* may be replaced 
by their average value of zero [33 such that probability 
amplitudes ai{t),a2{t) and a3{t) of the paired states are 
frozen approximately, effectiveness of the high-frequency 
approximation is lost partly. Particularly, in the case of 
moderate \u\ values, namely the values are neither very 
small nor too large, Eq. (4) is no longer a good approxi- 
mation. Anyhow, for a stronger reduced interaction with 
a larger \u\ value, we require to employ other approxima- 
tion methods and to explore the second order tunneling 
effects. 



III. SECOND-ORDER TUNNELING IN THE 

FAR-RESONANT STRONGLY-INTERACTING 

REGIME 

Now we consider the far-resonant case with a stronger 
reduced interaction to investigate tunneling dynamics of 
the system, by means of multiple-time-scale asymptotic 
analysis. In the high-frequency regime, we set e = J/uj 
as a small positive parameter and t' = ut is the rescaling 



Pi = 



E 



Jr'Az;) 
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P2 = 



Jn'(^)J-n'{^) 



n' — — oo 



!^+n' 



(8) 



for Uo/u + n' ^ 0. Therefore, Eqs. (6) and (7) are 
always definable and applicable except for the resonant 
case in which pi tends infinity. It is worth noting that for 
a stronger reduced interaction obeying |u| > J at least, 
the value of any term in the summations of Eq. (8) is 
less than e^^ such that e^pi may be a second-order quan- 
tity and Eqs. (6) and (7) could be applicable as a set of 
second-order equations. In fact, \pi\ < e~^ implies that 
the inequality \Uo/uj + n'\ = \n + n' + u/iu\ > e = J/lo 
holds for any pair {n G [0, cx)), n' G (—00, 00)}, which 
results in |u| > J. Particularly, we will numerically prove 
that perfect applicability of the second-order perturba- 
tion method requires |u| > lOJ later. Combining Eq. 
(6) with Eq. (7), we note that dynamics of the three 
paired states [the two bosons occupy the same site for 
Aj{t') with j = 1, 2, 3] is decoupled from that of the three 
unpaired states [the two bosons occupy distinct sites for 
Aj{t') with j = 4,5,6]. 

Clearly, for |u| > J, Eqs. (6) and (7) describe the sec- 
ond order approximation, where the time evolution of any 
paired state amplitude is a second order long-time-scale 
process, since its time derivative is proportional to only 
the second order constant e^. The similar results have 
been seen previously in the tight-binding optical lattice 



[26| . The second-order coupling coefficients of Eqs. (6) 
and (7) are proportional to the parameter e^/92, which 
describes the second-order tunneling rate of the system. 
The nonzero tunneling coefficient means that the tun- 
neling can occur, respectively, between the three paired 
states based on Eq. (6), and between the three unpaired 
states based on Eq. (7). In Fig. 2, we plot the factor p2 
of the second-order tunneling coefficients as a function 
of the driving parameters e/w and self-interaction inten- 
sity Uq/uj, where Fig. 2(b) is the plan view of Fig. 2(a). 
Combining Eq. (8) with Fig. 2 we can see that the fac- 
tor p2 tends to infinity for any integer value of Uq/uj and 
arbitrary value ranges of e/w, while its values are small 
enough for the considered far- resonant regime. Note that, 
in Fig. 2, the very great p values are not shown, since 
we have avoided the integer values of Uq/uj through se- 
lecting a rational step such that the multiple-time-scale 
asymptotic analysis holds. In the second approximation, 
Eqs. (6) and (7) mean that the CDT between the paired 
states will occur provided that the condition /92 = is 
satisfied. 



tions 




FIG. 2: (Color online) The factor p2 of the second-order tunneling 
coefficients as a function of e/oj and Uq/uj^ defined by Eq. (8), where 
Fig. 2{b) is the plan view of Fig. 2(a). The infinite p2 value at integer 
Uq/lo has been omitted. 
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Combining Eq. (9) with Eq. (11) produces ^^^ ^ = 

^^A^e"*^^*. Eliminating {A'^+A'^) from this equa- 
tion and Eq. (10) yields the decoupled equation 



zi; 



^:^i^ + .^A-o. (12) 



This is a second-order linear equation with constant co- 
efficients, whose general solution is well-known, A'2 = 
A'j^ exp(x+i) + A'_ exp(x_i) for the parameters x± = 

i{jjPi ± TTVPi+SpI) ^'^d the undetermined constants 
A'j. adjusted by the initial conditions. Under the above 
initial conditions the general solution becomes the special 
one 



Ait)^ 



cos(a;ii) — i 



Pi 



pI 



sin(wii) , (13) 



with LOi = —\/p\ + 8/92- Thus under the initial condi- 
tions ^2(0) = 1, Pj5^2(0) = 0, the analytical probabilities 
of Eqs. (6) and (7) are constructed as P4(i) — Psit) — 

P6it)=0, 



^2(i) = l^2(i)P = - 
Pi 



pI 



^pI 



p\ 



cos2(wit), (14) 



Pi{t) = Psit) = 21^0 2 ^in'{^it). (15) 

Secondly, for the two bosons initially occupying the left 



0.8 

0.6 


\t/ 


V: A 


a 0.4 


\'' J 


Y- '',/P4,5,6 


0.2 





AJ/. ■ 



50 100 150 200 




According to Eqs. (6) and (7), we make an exact com- 
parison of tunneling rates between the three paired states 
and three unpaired states for two different initial condi- 
tions as follows. Firstly, for the two bosons initially oc- 
cupying the middle well [i.e., P2(0) = 1, Pj^2{0) = 0], 
we seek the analytical solutions Pj{t){j — 1,2,..., 6) 
from Eqs. (6) and (7). To do so, we make the func- 
tion transformations Ai — ^' o^t^i'—i' ei 



Aiexp(-^^^i), A2 = 



A'2 exp(- 



■4JVl 



t) and A3 = A3 exp(- 



■2Jfpi 



FIG. 3: (Color online) Time evolutions of the probabilities Pj [j — 
1,2,..., 6) for the parameters J — 1^ uj — SO, e/tJ — 2, Uo/lj — 1.5 
and the initial conditions (a) ^2(0) ^ 1, Pj^2(0) ^ 0; (b) PsCO) ^ 
1, Pj^5(0) = 0. In (a), the dashed line corresponds to the probabilities 
Pi, 3, and the thin and the thick solid lines indicate the probabilities 
P2 and P4,5,6) respectively. In (b), the dashed line corresponds to the 
probabilities P^^q, and the thin and the thick solid lines indicate the 
probabilities P5 and Pi,2,3i respectively. In this figure and the following 
figures, all the circular points indicate the analytical solutions and the 
curves represent the numerical results. 



t). Inserting and middle well, respectively, [i.e., ^5(0) — 1, Pj^^{0) 



these expressions into Eq. (6) yields the coupled equa- 0], we easily obtain the analytical solutions Pi{t) 



P2{t) = P3(i) = 0, P5(i) = cos2(w2<), and P4(i) = 
PQ{t) = ^sin^(a;2i) with uj2 = V^JJoiz^), if '^^ neglect 
the second-order small quantities of Eq. (7) in the high- 
frequency regime. As an example, selecting the parame- 
ters w = 80 > J = 1, e/uj = 2 and Uq/lu = 1.5. There- 
fore, the tunneling period Ti — tt/wi « 89 corresponding 
to the second-order tunneling effect, and the tunneling 
period T2 = tt/lo2 ~ 10 corresponding to the first-order 
tunneling effect. 

In Fig. 3, we numerically plot the time evolutions of 
the probabilities Pj (j = 1,2, ...,6) based on Eq. (3) for 
the above two initial conditions, and the circular points 
correspond to the above analytical results. Obviously, 
the analytical results are in perfect agreement with the 
numerical simulations. The zero probability of the un- 
paired states in Fig. 3(a) and the zero probability of 
the paired states in Fig. 3(b) mean the selected CDT 
between the paired states and unpaired states. 



the same order of magnitude as the above Ti for a second 
order tunneling period. 
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FIG. 4: (Color online) Time evolutions of the probabilities Pj (j — 
1,2, ...,6) for the initial conditions Pi{0) = 1, Pj^4(0) = 0. The 
driving parameter is set as e/cj — 2.405, and the other parameters arc 
the same as those of Fig. 3(b). 

In the high-frequency regime, we note that the second- 
order correction term is much less than the first-order 
term, as shown in Eq. (7), and the dominating dynam- 
ics is decided by the first-order correction term provided 
that it does not vanish, on the contrary, is decided by 
the high-order correction term if the first-order correc- 
tion term vanishes in the case of Jo(^) = 0. As another 
example, in Fig. 4, we numerically show the time evo- 
lutions of the probabilities Pj [j = 1,2,..., 6) based on 
Eq. (3) for two bosons initially occupying the left and 
middle well, respectively [i.e., P4(0) — 1, Pj^a[Q) — 0], 
where the parameter is set as e/w = 2.405 correspond- 
ing to the first zero of J7o( — ) and the other param- 
eters are the same as those of Fig. 3(b). In this 
case, we readily calculate the probabilities analytically, 
because of A^{t) = in Eq. (7). Under the ini- 
tial conditions ^4(0) — 1 and Aj-t^^ifi) = 0, we im- 
mediately obtain the analytical solutions A\ — A2 = 
A3 = A5 = 0, yl4 = exp(i^^t)cos(2:^t) and A^ = 

exp(i — —t) sin( ^ ^^ t). Thus, the corresponding proba- 
bilities of Eqs. (6) "and (7) read P^ = P2 = P3 ^ P5 = 0, 

P4 = |A4(t)|2 = cosH^t) and Pe = sm^{Ml£it), 
which are plotted by the circular points in Fig. 4. The 
analytical and numerical results consistently verify that 
the tunneling period in Fig. 4 is about 130, which is in 
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FIG. 5: The time-averaged total probability (S) versus the self- 
interaction Uo in (a) and versus the driving frequency uj in (b), com- 
puted numerically based on Eq. (3). (a) the parameters are the same 
as those of Fig. 1 except for Uo, and (b) J — 1, e — 160, and Uo — 200. 
The time used for averaging is 200 in dimensionless units. 

In order to further confirm the analytical results in 
Eqs. (4), (6) and (7), we define the time-averaged to- 
tal probability of finding the two interacting bosons in 
the three paired states as (5) = (Pi) + (P2) + (P3) == 
7 /([(Pi + P2 + P3)dt for T = 200 J. The normalization 
means the time-averaged total probability in the three 
unpaired states being 1 — (S). Taking the initial con- 
ditions P2(0) = l,Pj^2(0) = and parameter J = 1, 
from Eq. (3) we numerically give (S) as the function 
of the self-interaction Uq for to = 80 and e/uj = 2.405, 
as in Fig. 5(a). It is shown that the time-averaged to- 
tal probability in the three paired states possesses dif- 
ferent features, for the multiphoton resonant points, the 
far-resonant regions and the near-resonant regions, re- 
spectively. Firstly, at each of the resonant points, i.e., 
Uq = muj with ?7i = 1, 2, ..., 5 being integer, {S) drops to 
the lowest points which mean that the separation proba- 
bility 1 — (S*) of the two bosons is the largest, as shown in 
Fig. 1. Secondly, the two bosons can also be separated in 
the near-resonant regions, however, the time-averaged to- 
tal probability {S) tends to one and the separation prob- 
ability tends to zero rapidly as increasing the reduced 
interaction strength \u\. Finally, in the far- resonant re- 
gions with larger |m| values, {S) is always equal to 1. Let 
half-width of the valley centred at mth resonant point 
of Fig. 5(a) be \u\m- The largest half-width for fitting 
{S) sa 1 can be estimated as \u\i k, lOJ from the first 
valley. This indicates that a selected CDT between the 
paired states and the unpaired states can happen in the 
region lOJ < |m| < a;/2, which is called the far-resonant 
strongly-interacting regime in this paper. Such a CDT 
enables the two bosons form a stable bound pair and 
cannot move independently for a stronger reduced inter- 



action. Similar to Refs. [l3|, ll8|, the phenomena can 
be understood that potential energy of two bosons oc- 
cupying a single well for strong repulsive interaction is 
greater than the maximum kinetic energy of two sepa- 
rate bosons, according to the principle of conservation of 
energy, the two bosons only forming a stable bound pair 
tunnel from a well to a neighboring well in the triple- 
well without dissipation. Only for the resonant case, the 
boson pair can be separated, because the bosons could 
absorb photons from the ac driving field. In the weakly- 
interacting regime, CDT is expected to occur in Fig. 5(a) 
because e/w = 2.405 is the first root of Jo{x) = 0, and 
we will further consider it in the next section. To show 
that the above analysis is generic in the high-frequency 
regime, we plot the time-averaged total probability (S) 
as a function of the driving frequency uj in Fig. 5(b) with 
£ = 160 and Uq = 200. Fig. 5(b) explicitly shows that 
the lowest points of (S) appear at the resonant points 
Ua/uj — 5,4,3,2,1 for the sufficiently high frequency, 
UJ > 30. The results agree with the above analysis on 
Fig. 5(a). 

It is well known that CDT can occur at the collapse 
points of the Floquet quasienergy spectrum [2J, |4l| , so 
the above-mentioned tunneling properties will be con- 
firmed by the Floquet quasienergy analysis as follows. 





FIG. 6: Numerical quasienergy spectrum versus e/uj for the self- 
interaction Uq — mcj, m — 0, 1, respectively. The parameters are 
set as J — 1, cj — 80, and (a) Uq — 0, (b) Uo — 80. The inset shows an 
enlargement of quasienergies near the point e/uJ — 2.405 corresponding 
to the first zero of ^^{e/uj). 



IV. FLOQUET QUASIENERGY ANALYSIS 

The Floquet theory provides a powerful tool to an- 
alyze the dynamics of a time-periodic quantum system 



|42l [. According to the Floquet theory, the solutions of 
the time-dependent Schrodinger equation can be written 
as IV'fc(t)) = e-'-^'=*|</.fc(t)), with \(l>k{t)) being the Flo- 
quet states and E^ Floquet quasienergies. In analogy 
to the Bloch solutions for the spatially periodic system, 
the quasienergy can only be determined up to a inte- 
ger multiple of the photon energy oj, and for the sake of 
definiteness it is usually assumed to vary in the first Bril- 
louin zone — a;/2 < E < ui/2. The Floquet states inherit 
the period of the Hamiltonian, and are eigenstates of the 
time evolution operator for one period of the driving 



U{T,0) =rexp 



H{t)dt 



(16) 



where T is the time-ordering operator and T = 2t: /lu is 
the period of the driving. Noticing that eigenvalues of 
U{T,0) are exp{—iEkT), the quasienergies of this sys- 
tem can be determined directly so long as we diagonalize 
C/(T, 0). In Fig. 6, selecting the parameters as J = 1 and 
ll) = 80, we show the numerical results of the quasienergy 
spectra as the functions of driving parameters s/uj for 
Uq = muj, 771 = 0, 1 with zero reduced interaction, re- 
spectively. In Fig. 6(a), for two noninteracting bosons, 
the quasienergy spectrum shows collapses at some fixed 
values of the driving parameters for which i7o(— ) = 0- 
The inset of Fig. 6(a) is an enlargement of quasienergies 
near the collapse point e/cu = 2.405, corresponding to the 
first zero of Jo{—) and shows an exact level-crossing at 
e/uj sa 2.405, analogous to a single boson in a triple-well 
system [29]. In the resonant regime with Uo = w = 80, 
the quasienergies of Fig. 6(b) show that the crossings of 
some quasienergies and the avoided crossing of the other 
quasienergies appear at the zero points (e/w = 3.832, ...) 
of the first-order Bessel function J'i{e/Ld). The numerical 
results are in good agreement with the analytical results 
from Eq. (4) with to = 1, u = in second section. In 
the case u = 0, Eqs. (6) and (7) are no longer valid and 
any quasienergy may be associated with both the paired 
states and the unpaired states. 

When the reduced interaction strength is sufficiently 
larger (e.g. |m| > J), the quasienergy spectrum is di- 
vided into two energy bands which correspond to the 
three paired states and the three unpaired states, respec- 
tively, as shown in Fig. 7, where the quasienergies of 
paired states (or unpaired states) aperiodically oscillate 
near u values (or value), so width of the energy gap be- 
tween the two bands is proportional to the |m| value. For 
a weaker interaction with Uq ^ u ^ 2, CDT between the 
different paired states and between the different unpaired 
states can be realized for the same driving parameters, as 
indicated by the level-crossing points in Fig. 7(a), when 
the ratio of the field amplitude e and the field frequency uj 
is a root of the equation Jo{e/u}) = 0. Precise agreements 
between the numerical results based on Eq. (3) and the 
analytical results from Eqs. (19), (20) and (23) are ob- 
served in Fig. 7(a) for a sufficiently larger range of ratio 
e/bj. The small deviation between the both results in Fig. 
7(a) indicates that the second-order perturbation method 
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FIG. 7: (Color online) Numerical quasicncrgy spectrum versus e/cj 
for (a.) Uo = u = 2J = 2 and (b) C/q = i^J - 30 = 50, u = -30. 
Tile other parameters are the same as tiiosc of Fig. 6. In Fig. 7(a), 
the above three curves denote the quasienergies of paired states and 
the below ones the quasienergies of unpaired states, while the situation 
is contrary to Fig. 7(b). The circular points are associated with the 
perturbation results from Eqs. (19), (20) and (23), the thin dotted 
lines indicate the zero points of J'o{e/uj) and the arrow in Fig. 7(b) 
indicates the amplification position of the inset. 



is perfectly applicable only for some suitable parameter 
regions. We have also investigated the quasienergy spec- 
tra for |u| = 5,6,... which are not shown here. The re- 
sults conformably verify that in the far-resonant strongly- 
interacting regime, uj/2 > \u\ > \u\i ~ lOJ is just the 
above suitable parameter regions. Interestingly, the en- 
ergy band corresponding to the paired states becomes 
narrower and the energy gap tends to wider as the in- 
crease of self- interaction intensity from \u\ = 2 < |u|i to 
|u| = 30 > |m|i. The wider gap means quantum transi- 
tion between the both kinds of states is hard to occur, 
and the narrower band necessitates to analyze the Flo- 
quet quasienergy spectrum from both cases of the un- 
paired states and the paired states, respectively. 



A. Avoided level-crossing of unpaired states 

In the far-resonant strongly-interacting regime, select- 
ing the parameters as J = 1, w = 80 and Uq = 50, 
i.e., m = 1, u ^ —30, from Eq. (3) we numerically 
plot quasienergy spectrum versus e/uj in Fig. 7(b). In 
this figure, the quasienergies corresponding to the un- 
paired states shows collapses when e/uj are the roots of 
Jo{e/L>j) = 0, however, the energy band corresponding 
to the paired states has collapsed into an approximate 
straight line. The inset of Fig. 7(b) is an enlargement of 
quasienergies corresponding to the three unpaired states 



near the first collapse point e/uj w 2.405, and the fine 
structure of energy spectrum exhibits that the pseudo- 
collapse point is converted to an avoided crossing point 
at e/oj K, 2.405 and two different crossing points due to 
the second-order correction terms in Eq. (7). 

To explain the numerical result, from Eq. (7) we ana- 
lytically calculate the quasienergies corresponding to the 
three unpaired states. Note that the period of func- 
tions exp[— i(^(i)] and exp[±2i(/3(t)] is T. Therefore, we 
can construct the Floquet states by setting |29i] Aj{t) = 
Bj exp{—iEt) [j = 4, 5, 6) for the three unpaired states 
with constant Bj, then rewriting Eq. (7) as the time- 
independent form 

£ ,P 

EBi = -JJ„{-)B5 - 2— (2B4P1 + Bep2), 

U! LU 

EB5 = -JJo{-){Bi + Be), 

U! 

£ 7^ 

EBe - -JJo{-)B5 - 2—{2BePi + B4P2). (17) 

U! LU 

The existence condition for the non-trivial solution of Eq. 
(17) reads 



E + A^p, JJo(f) 






-P2 



J Mi) E + i^p, 



(18) 



From Eq. (18) we obtain three Floquet quasienergies 
corresponding to the three unpaired states 



Ei = - 

E2 = - 

£'3 = — 
where we have set 



2(2JVi - J^p2) 

U! 
-2,Ppi - ,Pp2 - P3 

-2.Ppi - .Pp2 + P3 



P3 



{2J^pi + J^P2Y + 2.Poo^J§{-). 



(19) 



(20) 



We now compare the analytical results of Eqs. (19) and 
(20) with the numerical computation based on the origi- 
nal Eq. (3). A typical behavior of quasienegies near the 
first crossing point is plotted in the inset of Fig. 7(b). 
It is clearly shown that the analytical result (the circular 
points) is in perfect agreement with the direct numerical 
computation (the curves). 



B. A fine structure of quasienergy spectrum of 
paired states 

Next, we examine some detailed features of the 
quasiener gies corresponding to the three paired states. 
In Ref. [2^, Longhi et al proposed that in a lattice 
system, CDT can be realized between the paired states 



and between the unpaired states for the same parame- 
ters, namely the field parameters take the second root 
e/w = 5.52 of J'o{£/uj) — and the interaction inten- 
sity obeys Uq/uj — 2.58 corresponding to p2 — 0. Here 
for the triple- well system we prove the similar result, and 
exhibit a fine structure of quasienergy spectrum of paired 
states, based on analytical Floquet solutions of Eqs. (6), 
(7) and (8). According to Eq. (6), CDT occurs be- 
tween the paired states if the condition p2 = is sat- 
isfied. From Eq. (8), we have p2 = at e/uj « 0.95 
in the region e/uj e [0,8] for Uq/uj = 1.6, and have 
P2 = at e/uj w 1.20,2.02,5.52,5.74 in the same region 
for Uq/uj = 2.58. Selecting two different values of the 
self-interaction intensity, from Eq. (3) we numerically 
plot quasienergy spectrum versus e/uj in Figs. 8(a) and 
8(b), respectively, with the insets being enlargements of 
quasienergies of the three paired states. At the points 
fitting p2 — 0, the level-crossing of two quasienergies will 
occur, and this indicates that CDT for paired states can 
be observed at the crossing points of the partial levels. 
The predictions of the perturbation analysis can be con- 
firmed by direct numerical computation of the temporal 
evolution of the boson occupation probabilities from Eq. 
(3) (not depicted here). 



1, 2,3) into Eq. (6), we obtain easily 
.J' 



{E - Uo)Bi 



2—{Bipi+B2P2), 

UJ 
2 



,J 



{E ~ Uq)B2 = 2— [252/91 + (Bi + B3)P2] 
J2 

{E - Uo)B3 = 2 — (B3P1 + B2P2) 

UJ 



(21) 



with the existence condition of the non-trivial solution 



E-Uq 

J4 



4— pi 

UJ 



E-Uo 

J2 



J2 ,2 
2— Pi 

U! 



i—plfE-UQ-2—pA =0. 

W^ V U! / 



(22) 



From this equation we obtain three Floquet quasienergies 
corresponding to the three paired states as 



E^-- 


= Uq 

= Uo 

= Uq 


+ 2— Pi, 

UJ 

, 3J2pi + 






Kr 


VJ'p'i 


+ 8 JVi 


-f^5 


1 SJ'pi 


UJ 




i?fi = 


VJ'pi 


+ 8 JVi 



(23) 
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We note that the quasienergies Ej for j = 4,5,6 should 
be converted to E', in the first Brillouin zone, for example, 
Uq — 2.58a; w 206.4, so the quasienergies are rewritten 
as E'j = Ej - 3uj for w = 80 (j = 4,5,6). We plot the 
analytical quasienergies versus e/w for Uq/uj — 1.6 and 
2.58 in the insets of Figs. 8(a) and 8(b) as the circular 
points, respectively, which are in perfect agreement with 
the direct numerical computations (the curves) based on 
the original Eq. (3). 
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FIG. 8: (Color online) Numerical results of quasienergy spectrum 
versus e/u for the stronger reduced interaction case, (a) U^juj = 1.6, 
(b) Uq/u — 2.58. The other parameters are the same as those of 
Fig. 5. The arrows in each figure label the quasienergies of the paired 
states, and indicate the corresponding amplifications in the insets. The 
circular points in the insets are associated with the perturbation result 
from Eq. (23), and the thin dotted lines indicate the zero points of p2. 



Following, we analytically calculate the quasienergies 
corresponding to the three paired states. Substituting 
the Floquet solutions Aj{t) = Bj exp[— i(i? — C/q)^] (j = 



V. CONCLUSION 

We have investigated the tunneling dynamics of two 
bosons in a high-frequency driven triple well for a 
continuously-increasing interaction intensity, by means of 
the multiple-time-scale asymptotic analysis. In the ob- 
tained far-resonant strongly-interacting regime, we con- 
sider the second-order perturbed correction and show 
that the dominant tunneling effect of paired states is a 
second order process, similar to two bosons in a driven 
optical lattice [2^. For a stronger reduced interaction, 
we make an exact comparison of tunneling rates be- 
tween the paired states and unpaired states, and find 
that two bosons initially occupying the same well would 
form a stable bound pair. The selected CDT between 
the paired states and the unpaired states can occur for 
different values of interaction intensity. However, for the 
near-resonant case such initially paired bosons can sepa- 
rate due to the multiphoton resonance. Further we cal- 
culate the quasienergy spectrum and demonstrate that 
for the reduced interaction strength obeying \u\ > J, 



the quasienergy is divided into two energy bands corre- 
sponding to the three paired states and the three un- 
paired states, respectively. Width of the energy gap be- 
tween the two bands is proportional to the \u\ value. 
The prediction on the CDT is confirmed by the Floquet 
quasienergy spectra in which the avoided level-crossings 
and new level-crossings near the collapse points are ex- 
hibited for the three unpaired states, due to the second- 
order corrections. While for the three paired states, a fine 
structure of quasienergy spectrum up to the second order 
is displayed by which we show the different level-crossings 
beyond the former collapse points. The analytical results 
are very consistent with the direct numerical computa- 
tions from the time-dependent Bose-Hubbard Hamilto- 
nian for |m| > 10 J. The second order results of long time 
scale could be conveniently applied to adiabatic manipu- 
lationlsl, |43| of paired-particle tunneling in experiments 

Mm 
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Appendix: multiple-time-scale asymptotic analysis 

In the high-frequency regime, e = J/lj is a small pos- 
itive parameter. Let t' — ut be the rescaling dimension- 
less time variable, we rewrite Eq. (4) as 



.da\ 

.da2 

. das 
^ dt' 
.da^ 

^~dF 

'dt' 
,daQ 

''IF 



a4e[^^*'+^^(*')l-fa6e[^^*'-'^(*')l 



\/2e 

[-iv(t')i 

v(t')]" 



056 



a4e*'^(*') _,. agg[ 



+ 056 



ii^(t') 



(Al) 



with ip(t') = — sint'. At first, we transform the inde- 
pendent variable t' into the multiple-time-scale variables 
T„ = e"i', n = 0, 1, 2, ..., then replace the time derivatives 
by the expansion 



dV 



Oto + e^Ti + e^dr^ 



(A2) 



At the same time, we expand aj{t'){j — 1, 2, ..., 6) as the 
power series of e 

aj(i') = af\t') + eaf\t') + e\f\t') + ■■■. (A3) 

Substituting Eqs. (A2) and (A3) into Eq. (Al), and col- 
lecting the terms of the same order, we obtain a hierarchy 
of approximation equations of different orders in e. At 
leading order e°, one has 



da 



(0) 



dTo 



-0, af^=A,iT,,T2,---), 



(A4) 



where the amplitudes Aj(Ti,T2, ...) are functions of the 
slow time variables Ti, T2, ..., but are independent of the 
fast time variable Tq. At order e, one obtains the coupled 
equations 



.da\ 



(1) 



dTo 



.da. 



(1) 



. dAi 
'dT, 

dA2 



y2A4el'^^°-'^^(^«)l 



dTo 



dTi 

+ A(ie^'^'^°-"f^'^°'>^ 



V2 A4e[^^^°+'^(^'')1 



.da 



(1) 



dTo 



.da 



(1) 



. dAs 
'dTi 

.dAi 



y2A6e['^^°+*^(^°'l, 



V2Aie[-'^^°+'^(^")l 
dTo dTi V 

+ y2A2e[-^^^«-'^(^«)l -I- Ae[-'^(^°)1 



.di 



,(1) 



dTo 



.da, 



(1) 



dTo 



. dA^ 
'dTi 

.dA 
'dTi 



A4e^^(^«) + A6e[-'^(^°)1 



V2A2e^-'^^°+'^^^°'>^ 



+ V2A3e[-*^^°-^^(^°)l+Ae*^(^°) 



(A5) 



For the conveniences of our discussion, we simplify Eq. 



-(1) 



(1)/ 



(A5) as iday'/dTo = -idAj/dTi + Gy'{To) for j ^ 
1,2,..., 6. To avoid the occurrence of secular grow- 
ing terms in the solution a- , the solvability condition 
[24, 26] 






(A6) 



must be satisfied, where the overline denotes the time 
average with respect to the fast time variable To, i.e., 
the dc component of the driving term G, (To). The 
amplitudes a, at order e are given by 



4> ^ -i 



To r 



G^'\To) - G^'\To) 



d^. 



(A7) 
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Employing Eqs. (A5) and (A6), one gives 






■.dA4 



dTi " dTi 

dAa 



0, 



' dTi ~ " dTi 



;M5 






(A8) 



So the solutions of order e read 



a^"^ = V2iAiF„iTo), 

a'i'> = V2^[A4^l(To) + AeFoin)], 

a'i^ =V2iAeFi{To), 

ai^'> - i[V2AiF*{To) + V2A2F*{To) + A^i^a (^o)], 



^t[A4F2(To) + AeF*{To)] 

- i[V2A2F*{To) + V2A3F*{To) + A^F^i^ 



(A9) 



with 






FoiTo) 

FiiTo)=J2jn'{ 
F2iTo) 



^ CO -' 

'Uo 






y^ £ exp(mTo) - 1 



(AlO) 



n'i^ii 



We note that the probabilities to find the two strongly 
interacting bosons in the same wells are constants in 
time up to the first-order time scale T\ from Eq. (A8). 
Therefore, we need to consider the asymptotic analysis 
up to the order e^. Following the same procedure out- 



,(2 



(2) 



lined above, one has ida^' /dT^ ^ -~idAj/dT2+G)' (Tq), 



with 



(2) 



G\ 



G?) = 



V2iJo{-)A5Fo{To) - V2iiV2AiF*in) 

to 

+V2A2F*{To) + A^F^iTo)) cM^—To 

V2iM-)A^{Fa{To) + Fi{To)) - 
^/2^ [{V2AiF*{To) + A^F^in) + 

V2A2F*{To)) exp{i^To + i^{To)) + 
{V2A2F*iTo) + V2A3F*{T„) + 



Gf = ^/2^Jo(-)A5^l(To) - V2i(V2A2F* {T^) 
+V2A3F*{To) + AF2(ro)) eM^—To 

UJ 

+i(f{To)), 
Gf = 2Jo(-)(^4 + Ae)F*{To) - i\2A^Fo{To) 

UJ I 

exp(-i^To + i^{To)) + 2{A4F*iTo) 

CO 

+A6F*{To)) cxp{-i^To - lifiTo)) + 



G 



(2) 



{AiF2{n) + AeF*{To)) cxp{-iip{To)) 
iM-)A5{F2{To) + f;{To))- 

U! 

i[{V2AiF*{To) + V2A2F*{To) + 

A5F;iTo))exp{iipiTo)) + 
iV2A2F*{To) + V2A3F*{To) + 

A5^2(ro))cxpH^(To))l, 



G. 



(2) 



V2{A4Fi(To) 



iM-){A4 + Aq)F2{To) - t 

U! 

+A(iFn{n)) cxp{-i^To + MToj) + 

OJ 

(A4F2(ro) + AeF^iToj) cxp{i^{To)) ^ 
2A6^i(To)exp(-z^To - i^{n)) 

UJ 



(All) 



■^•^-«?' 



d 



'W2''' = ^^ 



(2) 



Then the solvability condition at order e^ reads 

'^ ^2{2Aipi + {Ai+Az)p2), 

'^ =2(A3Pl+A2P2), 

'^ =-2(2^4/91+^6^2), 

^ = 0, 

2T 



i-TT^A:i = G3 



d_ 
._d_ 
._d_ 
. d 



(2) 



l-TT^Ai = G4 



(2) 



l-TT^A^ = G5 



(2) 



.(2) 



^--A6 = G\^' = -2{2Aqpi + A4P2), 
0I2 

where we have set 



, ^ +n' 



Uo 



(A12) 



(A13) 



for Uq/lu + n' 7^ 0. Thus the evolution of the amplitudes 
Aj up to the second-order long time is given by 



dAj _ d 



-A^^'W2^^^' ^^''^ 



from Eq. (A3) and Eq. (A4). The corresponding proba- 
bility amplitudes read 



aAt')^A,{t')+o{e)+o{e^) + 



(A15) 
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in the high-frequency regime, where the high order small 
terms can be neglected. Substituting Eqs. (A4), (A8) 



and (A12) into Eq.(A14), we obtain the two sets of cou- 
pled equations, Eqs. (6) and (7). 
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